
import numpy as np
import matplotlib.pyplot as plt

# 定义参数
r1 = 4.95  # 圆柱半径
rc = 0.5  # 圆柱半径
t = 0.1  # 壁厚
E = 1  # 弹性模量 (单位: Pa)
mu = 0.3  # 泊松比

s0 = 1
s1 = 0
a=0.1

# 定义函数
def uy(s0, rc, fai, E):
    return s0 * rc * np.sin(fai) / E

def ur(s1, v, rc, fai, E,mu_):
    return s1 * (1 + v) * rc * np.sin(fai) / E-0*s0 * rc * np.sin(fai) / E


def X(rc, fai):
    return rc * np.cos(fai)

def Y(rc, fai):
    return rc * np.sin(fai)

def Z(fai, rc, rl):
    return -rl + np.sqrt(rl**2 - rc**2 * np.cos(fai)**2)

def u_x_2D(ur_, uy_,rl, X_,Z_,mu_):
    return ur_
    return ur_ * (rl+Z_) / np.sqrt((rl+Z_) **2 + X_**2)

def u_y_2D(ur_, uy_, rl,X_, Z_,mu_):
    return uy_

def un(ux2d, uy2d, fai):
    #return uy2d * np.sin(fai)
    return ux2d * np.cos(fai) + uy2d * np.sin(fai)

def ut(ux2d, uy2d, fai):
    #return uy2d * np.cos(fai)
    return -ux2d * np.sin(fai) + uy2d * np.cos(fai)

# 参数化角度
fai = np.linspace(0, 2 * np.pi, 200)
fai=fai[:-1]

# 计算 X, Y, Z
X_ = X(rc, fai)
Y_ = Y(rc, fai)
Z_ = Z(fai, rc, r1)

# 计算 uy 和 ur
uy_ = uy(s0, rc, fai, E)
ur_ = ur(s1, mu, rc, fai, E,mu)

# 计算 u_x_2D 和 u_y_2D
u_x_2D_ = u_x_2D(ur_, uy_,r1,X_, Z_,mu)
u_y_2D_ = u_y_2D(ur_, uy_,r1, X_, Z_,mu)

# 计算 un 和 ut
un_ = un(u_x_2D_, u_y_2D_, fai)
ut_ = ut(u_x_2D_, u_y_2D_, fai)

print(un_)
print(ut_)

# 创建图形
plt.figure(figsize=(12, 6))

# 绘制 un 与 fai 的关系
plt.subplot(1, 2, 1)
plt.plot(fai, un_, label='$u_n(\phi)$', color='blue')
plt.xlabel(r'$\phi$ (radians)')
plt.ylabel(r'$u_n$ (m)')
plt.title(r'Relationship between $u_n$ and $\phi$')
plt.legend()
plt.grid(True)

# 绘制 ut 与 fai 的关系
plt.subplot(1, 2, 2)
plt.plot(fai, ut_, label=r'$u_t(\phi)$', color='red')
plt.xlabel(r'$\phi$ (radians)')
plt.ylabel(r'$u_t$ (m)')
plt.title(r'Relationship between $u_t$ and $\phi$')
plt.legend()
plt.grid(True)

# 显示图形
plt.tight_layout()
plt.show()

import k
print(k.predict_xy_edge_crack_2d(E,rc*2,a,un_,ut_,8))

import numpy as np

def calculate_stress_intensity_factor(sigma_2, D, t, a):
    """
    计算薄壁管道周向裂纹的应力强度因子 K_I。

    参数:
    p : float
        管道内压 (单位: Pa)
    D : float
        管道外径 (单位: m)
    t : float
        管道壁厚 (单位: m)
    a : float
        裂纹半长 (单位: m)

    返回:
    K_I : float
        应力强度因子 (单位: Pa·m^(1/2))
    """ 
    # 计算鼓胀因子 M
    M = np.sqrt(1 + 0.64 * (a / (D * t)))
    
    # 计算应力强度因子 K_I
    K_I = M * sigma_2 * np.sqrt(np.pi * a)
    
    return K_I

# 示例使用
D = r1*2  # 外径 10 cm
K_I = calculate_stress_intensity_factor(s0, D, t, a)
print(f"应力强度因子 K_I: {K_I:.4f} Pa·m^(1/2)")